function example4

% Solves the nonlinear BVP
%	y'' = f(x, y, y')  for  xL < x < xr  '
% where
%	y(xl) = yL  and  y(xR) = yR

% clear all previous variables and plots
clear *
clf

% set boundary conditions
xL = 0.0;	yL = 0.0;
xR = 1.0;	yR = 0.0;

% parameters for calculation
n = 6;
tol = 0.000001;

% calculate solution
x=linspace(xL,xR,n+2);
y=newton(x,yL,yR,n,tol);

% plot computed solution
subplot(2,1,1)
plot(x,y,'--ro','LineWidth',1,'MarkerSize',7)
hold on

% define axes used in plot
xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
title(['Nonlinear BVP: \mu = 1'],'Color','k','FontSize',14,'FontWeight','bold');

% plot exact solution
cc=0.7585822995;
xe=linspace(xL,xR,45);
for ix=1:45
	exact(ix)=log(2*(cc^2)*(1-tanh(cc*(xe(ix)-1/2))^2));
end;
plot(xe,exact,'k','LineWidth',1)
legend([' N = ', num2str(n)],' Exact','Location','South')

% have MATLAB use certain plot options (all are optional)
box on
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold'); 
axis([0 1 0 0.1501]);
hold off


% iteration to determine the error as number of points used is increased
ii=0;
n=1;
tol = 0.000000001;
for i=0:4
	n=5*n;
	ii=ii+1; points(ii)=n;

	% calculate solution
	x=linspace(xL,xR,n+2);
	y=newton(x,yL,yR,n,tol);

	for ix=1:n+2
		exact2(ix)=log(2*(cc^2)*(1-tanh(cc*(x(ix)-1/2))^2));
	end;
	
	% calculate error
	errorM(ii)=norm(exact2-y,inf);
end;

% plot error curve
subplot(2,1,2)
loglog(points,errorM,'--or','LineWidth',1,'MarkerSize',7)
grid on

% define axes used in plot
xlabel('Number of Grid Points','FontSize',14,'FontWeight','bold')
ylabel('Error','FontSize',14,'FontWeight','bold')

% have MATLAB use certain plot options (all are optional)
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14);
axis([1 1e4 1e-9 1e-3]);
set(gca,'ytick',[1e-9  1e-7 1e-5 1e-3]); 
hold off


function g=f(x,y,z)
g=-exp(y);

function g=fy(x,y,z)
g=-exp(y);

function g=fz(x,y,z)
g=0;


% Newtons method for solving the nonlinear system  
function y = newton(x,yL,yR,n,error)

% start off with a linear solution
% y=linspace(yL,yR,n+2);

% start off with a quadratic solution
mu=1;
for ix=1:n+2
	y(ix)=-mu*x(ix)*(1-x(ix));
end;

dx=x(2)-x(1);
dxx = dx*dx;
err=1;
counter=0;
while err > error

	% calculate the coefficients of finite difference equation
	a=zeros(1,n); c=zeros(1,n); v=zeros(1,n); u=zeros(1,n);
	for j = 2:n+1
		jj=j-1;
		z = (y(j+1) - y(j-1))/(2*dx);
		a(jj) = 2 + dxx*fy(x(j), y(j), z);
		c(jj) = -1 - 0.5*dx*fz(x(j), y(j), z);
		v(jj) = - 2*y(j) + y(j+1) + y(j-1) - dxx*f(x(j), y(j), z);
	end;

	% Newton iteration
	v(1) = v(1)/a(1);
	u(1) = - (2 + c(1))/a(1);
	for j = 2:n
		xl = a(j) - c(j)*u(j-1);
		v(j) = (v(j) - c(j)*v(j-1))/xl;
		u(j) = - (2 + c(j))/xl;
	end;
	vv = v(n);
	y(n+1) = y(n+1) + vv;
	err = abs(vv);
	for jj = n:-1:2
		vv = v(jj-1) - u(jj-1)*vv;
		err = max(err, abs(vv));
		y(jj) = y(jj) + vv;
	end;
	counter=counter+1;
end;
%newton_iterations=counter